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ABSTRACT 

Using high-resolution, two-dimensional hydrodynamic simulations, we investigate 
nonlinear gravitational responses of gas to, and the resulting drag force on, a very 
massive perturber Mp moving at velocity Vp through a uniform gaseous medium of 
adiabatic sound speed a^o. We model the perturber as a Plummer potential with 
softening radius r^, and run various models with differing A = GMp/ {a^rg) and A4 = 
Vp/ooQ by imposing cylindrical symmetry with respect to the line of perturber motion. 
For supersonic cases, a massive perturber quickly develops nonlinear flows that produce 
a detached bow shock and a vortex ring, which is unlike in the linear cases where 
Mach cones are bounded by low-amplitude Mach waves. The flows behind the shock 
are initially non-steady, displaying quasi-periodic, overstable oscillations of the vortex 
ring and the shock. The vortex ring is eventually shed downstream and the flows 
evolve toward a quasi-steady state where the density wake near the perturber is in 
near hydrostatic equilibrium. We find that the detached shock distance 5 and the 
nonlinear drag force F depend solely on 77 = A/{M'^ — 1) such that 5/rs = rj and 
F/Fiin = {rj/2)~^''^^ for rj > 2, where Fn^ is the linear drag force of Ostriker (1999). The 
reduction of F compared with Fn^ is caused by front-back symmetry in the nonlinear 
density wakes. In subsonic cases, the flows without involving a shock do not readily 
reach a steady state. Nevertheless, the subsonic density wake near a perturber is close 
to being hydrostatic, resulting in the drag force similar to the linear case. Our results 
suggest that dynamical friction of a very massive object as in a merger of black holes 
near a galaxy center will take considerably longer than the linear prediction. 

Subject headings: black hole physics — hydrodynamics — ISM: general — shock waves 



INTRODUCTION 



A massive object in orbital motion suffers from orbital decay due to a negative torque caused 
by gravitational interaction with its own gravitationally induced wake created in the background 
medium. This process, com monly referred to as d ynamical friction (DF), occurs in not only col- 
lisionless environments (e.g., IChandrasekhaiill943l ) but also collisional gaseous backgrounds (e.g., 
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Ostrikeiill999l ) . The DF in a gaseous medium is of great importance in understanding the formation 
and evolution of planets, binary stars, supermassive black holes, etc. For instance, gravitational 
interaction between a protoplanet and its environmental disk causes the former to migrate toward a 
central star, natura lly explaining the presence of "hot Jupiters" found from radial velocity surveys 
([Butler et al.l 120061 . and references therein). The migration also helps a planet grow faster in mass 
by providing an expanded gas-feeding zone at an enhanced accretion rate, which may overcome the 
failu re of in situ core-ac cretion scenario in building a giant planet within a typical disk lifetime (see, 
e.g., lAlibert et al.ll2005l ). In the case of nuclear black holes in merging galaxies, they are expected 



to first experience the DF to form a binary and then coalesce into a supermassive black hole by 
emitting gravitational waves. Friction of nuclear black holes against the collisionless stellar back- 
ground appears to be inefficient due to scatt ering and depletion of stars ne ar the black holes, which 
is known as the "final-parsec problem" (see iMilosavljevic MerrittI 120031 . and references therein). 
However, recent numerical A^-body/SPH simulations show that the gravitational drag from the 
gaseous ba ckground is sufficient to form a black-hole binary in a relatively short time (^ 1 — 10 



Myrs) fe.g.. lEscala et al.ll2004l . l200,4lDotti et al.ll200fil . l2007l : lMaver et al.ll2007l : ICuadra et al.ll2009l ) 



Thanks to a seminal paper of lOstrikeii (|l999l ). DF in a gaseous medium is well understoo d 
as long as density w a kes h ave si nall amplitudes. Ear l ier th eoretical work by iDokuchaevI (| 19641 ) , 
Ruderman Sz Spiegell (|l97ll ) , and iRephaeli Salpeteii (jl98d ) considered density wakes in a steady 
state and found that the drag force vanishes for a subsonic perturber, while it becomes remark- 
ably similar to the collisionless counter part for s upers onic cases. Using a time-dependent linear 
perturbation theory, on the other hand, lOstrikeri (| 19991 ) found that the gravitational drag force on 
a point-mass perturber with mass Mp moving at velocity Vp on a straight-line trajectory through 
a uniform gaseous medium with density poo and sound speed Ooo is given by 

ATTp^jGMpf j llni^)-M, M<1, 
'2Hl-j^)+H^), M>1, 



y2 



(1) 



where A4 = Vp/aoo is the Mach number, t is the time elapsed after the introduction of the perturber, 
and Tmin is the minimum radius introduced to avoid the singularit y in the force evaluatio n. Equa- 
tion ([1]) shows that the gaseous DF force becomes identical to the IChandrasekhaii (| 19431 ) formula 
for the collisionless drag for Ai ^ 1, and is, albeit small, non-zero even for su bsonic perturbers be - 
cause the time dependency breaks the symmetry in the density wakes (see also I Just &: Kegellll990l ). 
Equation ([T|) has been applied t o various astrophysical situations inc luding orbital decay of compact 

and heating of an intracluster 



objects in accretion disks (e.g., iNaravanI l2000l : iKaras &: Subi 



medium by supersonically moving galaxies in clusters (e.g. 



2001 



El- Z ant et al. 



2005; Kim et al. 2005 



Kim 



20071 : IConrov k Ostrikerll2008l l 



20041 : 



Faltenbacher et al 



While the result of lOstrikeii (|l999l l is valid in a strict sense only for a linear-trajecto ry perturber 
i n a un iform medium, it has proven to be applicable to more general cases. For example. iKim &: Kim 
(j2007l ) considered a circular-orbit perturber with orbital radius Vp in a uniform gaseous medium 
and found t hat equation (HI) is a reasonable approxim ation to the gaseous drag force on it, provided 



Vpt 



2r. 



Sanchez-Salcedo &: Brandenburg (j200ll ) numerically found that the orbital decay of 
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a Plummer sphere with a softening radius in a radially-stratifi ed medium is co nsistent with 
the prediction of equation ([1]), if Vpt/r^ia = (0.35rp/rs)^-^'*. Also, iBaraussd (j2007l ) showed that 
equation ([T|) remains valid even for a perturb e r with relativistic speed if the relativistic correction 
factors are included. While iDotti et al.l (|2006l ) found that the orbital decay of black hole binaries 
took longer than the prediction of equation ([T|) for a single perturber, the discrepancy between 
the numerical and analytical results can be reconciled, at least partly, by taking allowance for fact 
that an object in a binary experiences not only a negative torque d ue to its ow i i wak e but also a 
positive torque from the companion wake. For an equal-mass binary, iKim et al.l (120081 ) found that 
the positive torque is on average about 40% of the negative torque. 



Since the results of lOstrikeij are based on the assumption that density wakes remain in the linear 
regime, the validity of equation ([T]) for very massive perturbers has yet to be seen. The strength 
of gravitational perturbations due to a body with mass Mp can be measured by the dimensionless 
parameter 



A 



GMr, 



(2) 



which roughly corresponds to the perturbed density at a distance from the perturber relative 



to the background density (e.g., lJust Kegell Il990l : lOstrikeii Il999l ). and is equal to the Bondi 
radius r-Q = GMp/a^ relative to r^. For systems with A ^ 1, the density wakes are clearly in 
the nonlinear regime and the linear perturbation analyses are likely to fail. Identifying with 
the gravitational softening radius of a perturber (or, equivalently, its size), A is in the ranges of 
~ 0.1 — 1 for galaxies embedded in typical intracluster media, ~ 10 — 100 for protoplanets in 
protostellar disks, and ~ 10^ — 10^ for supermassive black holes near galaxy centers, sug g esting 
that the wakes of massive compact objects can readily be nonlinear. Indeed, lEscala et al.l (120051 ) 
reported that the orbital decay time of supermassive black hole binaries with ^ ~ 10 — 100 depends 
on Mp much less sensitively than the results of the linear theory, which may be caused primarily 
by the nonlinear effects. 

In this paper, we investigate nonlinear DF of a very massive perturber in a gaseous medium 
using numerical hydrodynamic simulations. In order to isolate the effects of the perturber mass and 
its velocity on the DF force, we c onsider a p ertur ber following a straight-line trajectory in a uniform 
gaseous medium, similarly to in lOstrikerl (119991 ). We model the perturber as a Plummer sphere 
that does not possess any solid surface and merely provides gravitational potential perturbations 
to the background medium that would otherwise remain static and uniform; to make contact with 
the results of the linear theory, we ignore the accretion of gas onto the perturber in the current 
work. Our primary objectives are to find the changes in distributions of density wakes with A and 
Ai, and to quantify the resulting gravitational drag forces in comparison with the linear cases. 

Nonlinear responses of a background to a massive perturber moving at a superso nic speed have 
been extensively studied in the context of the Bondi-Hoyl e-Lyttleton (BHL) accretion (iHoyle &: Lyttleton 
I939I : bondi fc Ho^ll944l : lBondll952l : see also review of lEdea and references therein) . Unlike 
our models where mass accretion to a perturber is prohibited, however, the BHL accretion problem 
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consid ered a perturber c ontaining a defined surface through which gas is either accreted or reflected. 
It was iHuntI (|l97ll . Il979l ) who first solved for BHL accretion fiows numerically, finding that the col- 
lisional nature of gas supports a bow shock in front of a supersonic perturber. Later studies found 
that the BHL accretion fiows exhibit unstable behaviors such as fiip-flop mo tions of the accretion 



shocks and vortex shedding when the condition of axisymmetry is relaxed (e.g.lMatsuda et al 



Frvxe 
1997I . 



1 &; Taam 
I999I : 



1988- 



Foglizzo et al. 



Taam &: Frvxell 

2nnd ). 



1988 



Matsuda et al. 



19891 : 



Ruffert 1994 



1987 



Foglizzo fc Ruffert 



While these numerical works on the BHL accretion explored 
temporal evolution and distribution of density wakes as well as nonlinear features in some great 
detail, they mainly concentrated on the gravitatio nal focusing and resulting accretion rate of gas 



onto the perturber. Although some authors (e.g., IShima et al 



1985 



Shankar et al 



1993 



Ruffert 



19961 ) presented values for aerodynamic and gravitational drag forces, only a limited range of A was 
covered. In this work, we run a number of numerical simulations by varying A and A4 systemati- 
cally in order to quantify the dependences of the gravitational drag force on these parameters. A 
brief comparison between our results with those from the BHL accretion studies will be presented. 

This paper is organized as follow: In ^ we describe numerical methods we employ for nonlinear 
simulations. In J}3l as a code test we revisit the cases with a spatially-extended, linear perturber 
with ^ = 0.01. We compare the resulting distributions of density and velocity wakes with those of 
analytical results and provide a way to handle the effect on the DF force of a softening radius which 
is necessary for numerical simulations. Evolution and quasi-steady distributions of fully nonlinear 
density wakes and the associated drag forces are presented in ^ Finally, in ^we summarize our 
findings and discuss their astrophysical implications. 



2. NUMERICAL METHOD 

In this paper, we consider an unmagnetized, inviscid, non-radiating, non-self-gravitating gaseous 
medium and study its gravitational responses to a massive perturber moving along a straight-line 
trajectory using numerical simulations. The background gaseous medium is initially static and uni- 
form with density poo and adiabatic sound speed Ooo • We adopt an adiabatic equation of state with 
an index 7 = 5/3 throughout simulations. The simulations are carried out on a two-dimensional 
(-R, z) plane in cylindrical symmetry, where R and z denote the distances from and along the axis 
of symmetry, respectively. We assume that the perturber exerts only gravity to the surrounding 
medium and does not possess any surface, so that neither accretion nor reflection of the gas is al- 
lowed. The perturber is modeled as a Plummer sphere with mass Mp that is moving at a constant 
speed Vp along the R = axis toward the positive z-direction, with the gravitational potential 

where is the softening radius. 

We take r,, clqq, and tcross — '"s/'^oo as the units of length, velocity, and time, respectively, in 
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our simulations. Then, our models are completely parametrized by A and M. We run a total of 
58 models with A varying from 0.01 to 600 and M in the range of .5 to 4.0. We solve the basic 
equations of ideal hydrodynamics using FLASH3 ([Frvxell et al.ll2000l ) . an Eulerian hydrody i iamic s 
code that implements a direct Piecewise-Parabolic Method solver of IColella &: Woodward! (|1984| ) 
for high-order spatial reconstruction. Although the FLASH3 code is capable of both uniform grid 
and adaptive mesh refinement calculations, we adopt the uniform grid method since the accurate 
evaluation of the drag force requires the whole computational domain to be well resolved. As we 
will show below, we find that it is necessary to have at least 5 zones per to obtain converged 
results for the drag forces. Our largest grid models have 3,072x12,288 zones in (R,z); we make 
sure that our computational domain is taken to be large enough to contain the whole density wake 
in a given model. The simulations are typically carried out until t/tcross = 600 when most of the 
wakes are well resolved and reach a quasi-steady state. 



3. LINEAR CASES 



Time-dependent linear perturbation theories for the DF drag force in a gaseous medium usually 
study the responses of gas to a low-inass, point-mass pe rturber corresponding to = in the 
Plummer potential (e.g.. Just &: Kegei 1990 : Ostriker 1999 ). which requires to introduce the cut-off 
radius rmin in the linear force formula (eq. [1]). In numerical simulations, on the other hand, one 
needs to assign a non-zero value to the softening radius, which in turn makes it unnecessary to 
use the cut-off radius in the force evaluation. Since our goal is to compare the nonlinear drag 
force on a massive perturber with the linear prediction, we have to first find a proper relationship 
between and r^[^ that makes the numerical and analytical results consistent with each other when 
j4 <C 1. Motivated by this consideration, in this section we briefly present the results of numerical 
simulations for a low-mass perturber with A = 0.01, and compare the resulting distributions of 
density and velocity wakes and the drag forces with those from the linear theories. This will also 
allow us to check the accuracy of our numerical experiments. 

Regarding the linear wakes with which numerical results will be compared, it is worth men- 
tioning that th ere are several analy tical methods for finding solutions for the perturbed density and 
velocity fields. lJust &: Kegell (|l99d ) utilized Fourier transform for the space variables and Laplace 
transf orm for t h e tini e variable, finding expressions both for the density and velocity wakes. In- 
stead, lOstrikeii (119991 ) used a retarded Green's functi on technique a nd fo und an expression only 
for the density wake that is identical to the result of Just &: Kegel (jl99ol ). We found that while 
the analytical formula for the density wake agrees well with our numerical results for a low-mas s 
perturber, the simulated velocity field differs from the expression given by lJust &: Kegell (jl990l ). 
In Appendix [XI we revisit the time-dependent linear theory using Fourier transforms both for the 
space and time variables. As we will show below, our expressions ()A19P and ()A20P for the per- 
tu rbed velocities are in good agreement with the numerical results, confirming that equation (47) 
of I Just Sz Kegell (jl990l ) contains a typographical mistake. 
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Figure [T] shows as color-scale images the snapshots on the z-R plane of the perturbed density 
a = pI Poo — 1 (top)) the parallel velocity (middle), and the perpendicular velocity vr (bottom) 
to the line of motion for a model with A = 0.01 and Ai = 1.5. Note that the coordinates are 
normalized by aoot- The perturber initially introduced at (z, R) = (0, 0) has moved to (Maoot, 0) 
at time t. The characteristic features of a supersonic wake consisting of a sonic sphere with radius 
a^ot centered at the initial perturber location and a Mach cone bounded by Mach waves located at 
R = —(Ai"^ — l)~^/^(z — Aittoot) for z > aoot/M. are apparent in the top panel. Also plotted as 
black solid contours are the results of the linear perturbation theory (eqs. [AlOj . |A19j . and |A20j ). 
which are overall in good agreement with the numerical results. A careful inspection of Figure [TJ 
however, reveals that the numerical results deviate slightly from the analytical ones especially near 
the sonic radius and the Mach waves. This is due to the fact that the perturber in our numerical 
models is modeled as an extended Plummer sphere rather than a point mass. 

One can semi-analytically construct the density wake aext of an extended perturber by convolv- 
ing the density wake a d ue to the corresponding point mass with the extended mass distribution 



Pcxt of the perturber (e.g.. lJust &: Kegellll990l : lFurlanetto &: Loebll2002l ). For a Plummer sphere we 



use, the convolution theorem gives 

acxt(x,t) = j a(x - x',t)pext(x',t) dV, (4) 

where pe:^t(x, t) = 2>Mprl(R^ + rl + (z- Vpt)"^)-^/"^ / (An). Figure plots as solid lines the profiles 
of CKext for A = 0.01 and Ai = 1.5 along the cuts at R/a^ot = 0.20 and z/aoot = 0.92 marked 
as dotted lines in the top panel of Figure [H which are in excellent agreement with the simulation 
outcomes (open circles). Compared with the point-mass results (dashed lines), the extended mass 
distribution tends to smear out the discontinuities at the boundary of the sonic sphere and the 
Mach cone. This makes sense since a perturbed density at one location is a superposition of sonic 
perturbations with various strengths launched by all the mass elements comprising the extended 
body. 

Allowing for the extended mass distribution, the gravitational drag force exerted on the per- 
turber in the negative z-direction can be obtained by directly evaluating the integral 

F(t) = II gPext(^VM^,t)(.-/) ^3^,^ 

where p(x, t) is the wake distribution. Although equation ([5]) requires the ai'-integration to be 
performed over the entire Plummer sphere, we empirically found that the drag force on the region 
with distance from the perturber larger than lOr^ has a negligible contribution to the total. Thus, 
in practice, we limit the integration to the region within lOr^ that contains about 98.5% of the total 
perturber mass. Figure [3] plots as open squares the numerical drag forces on a low-mass perturber 
with A = 0.01 but differing Ai. The solid line corresponds to equation ([T|) with 

r^in = 0.35>[0-V„ (6) 
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which gives the best fit to the supersonic results of our adiabatic simulations. Equation ^ is our 
prescription for t he cut-off radius when we compare the numerical drag forces with the analytic 
results. Note that ISanchez-Salcedo &: Brandenburd (jl999l ) suggested rmin = 2.25rs based on their 
isothermal simulations. The difference between the two prescriptions may be due in part to using 
differe nt equations of state and in part to low resolution (1 cell per r^) in lSanchez-Salcedo &: Brandenburg 



4. NONLINEAR CASES 
4.1. Wake Evolution 

4.1.1. Supersonic Cases 

We begin by describing the temporal evolution of our fiducial model with A = 20 and A4 = 
1.5; the evolution of other supersonic models are qualitatively similar. Figure H] illustrates the 
density and velocity structures of this model in a comoving frame with the perturber located at 
{s = z — Vpt,R) = (0,0). Only the region with \s\/rs < 20 and < R/rg < 20 is shown. Unlike 
the linear cases with ^ ^ 1 where the density wakes have too small amplitudes to launch shock 
waves, the perturber with ^ = 20 emits strong perturbations that quickly develop into a bow 
shock. The upstream gas moving toward the perturber along the symmetry axis is first accelerated 
by the gravity of the perturber, is shocked to subsonic speed, and then piles up near the center 
of the perturber. This creates a steep pressure gradient in between the perturber and the shock, 
tending to push the shock away from the perturber. On the other hand, the gas flowing above 
(not far away from) the symmetry axis is deflected toward the perturber even before entering the 
shock and decreases its speed after the shock. This gas thus has a longer time to be exposed to 
the gravity of the perturber as it moves toward the symmetry axis. The gravitational potential 
is so deep that the material arriving at the rear side of the perturber can be pull ed back toward 



the perturber, creating a stagnation point just as in the BHL accretion flows (e.g., iMatsuda et al 



19871 : iFryxell et al.l 119871 ) . This produces a strong counterstream that moves into the upstream 
direction along the symmetry axis, as well as a primary vortex in the s~R plane (Fig. IDa jl]. The 
counterstream combined with the pressure gradient in the front side pushes the shock front away 
from the perturber. 

The advance of the shock front in the upstream direction allows more time for the shocked gas 
to be affected by the gravity, strengthening the counterstream (as well as the primary vortex) and 
thus increasing the detached distance 6 of the shock measured along the symmetry axis. At the 
same time, the center of the vortex with a lower density than the surrounding buoyantly rises toward 



^More precisely, the vortex in the s-R plane is a vortex ring in three dimensions, and the counterstream is a part 
of the vortex ring near the symmetry axis. 
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the high-i? regions. This decreases the ram pressure of the counterstream exerted on the shock 
front and the shock advance slows down. Figure [5] plots the time evolution of 5 for some selected 
models, while Figure [6] traces the trajectory of the center of the primary vortex in our standard 
model [A = 20 and M = 1.5) on the s-R plane. For our standard model, 5 keeps increasing to ISr^ 
at t/tcross = 100. At this time, the vortex is located near at {s,R) ~ (— 0.2rs,5rs) (Figs. [Da and 
[6]). The shock front soon overshoots a potential equilibrium position where the thermal pressure 
of the postshock gas supports it against the perturber gravity, and begins oscillating around the 
equilibrium position. 

The shock oscillation changes the velocity field in the postshock region, causing the vortex 
to move in the counterclockwise direction around its mean position (— 0.5rs,5rs) following the 
background flow. Kelvin's circulation theorem implies that the vortex becomes stronger as it 
moves toward the symmetry axis, amplifying the strength of the counterstream near the perturber. 
The velocity of the counterstream becomes largest when the vortex arrives closest to the perturber 
on its oscillation path (Fig. HH). The shock front that was moving toward the perturber is pushed 
by the strong counterstream, and reverses its motion. The vortex rises to high R as the shock 
moves away from the perturber, weakening the counterstream, and the oscillation cycle repeats 
quasi-p er io dically. 

Figures [5] and [6] show that the amplitudes of the detached shock oscillations as well as the vor- 
tex movements grow secularly with time as the shock oscillation continues. This can be understood 
as follows. When the shock front is displaced from the equilibrium position away from the per- 
turber (e.g.. Fig. Hlc), the shocked subsonic gas can acquire an extended time to be gravitationally 
influenced, similarly to the early situation when the shock advances away from the perturber. With 
the stagnation point moving away from the perturber, more mass and momentum are added to 
the counterstream, amplifying the vortex oscillation. In addition, when the strong counterstream 
collides directly with the shocked subsonic flow near the symmetry axis, the material at the inter- 
face is injected upward in the lateral direction and then carried in the negative s-direction by the 
background flow (see, e.g., the region at s/r^ ~ 12 and R/rg ~ — 3 in Fig. [4]e). This produces 
small vortices near the interface that move downstream and merge with the primary vortex, again 
strengthening the latter. Consequently, the primary vortex is able to move closer to the symmetry 
axis and amplify the counterstream in the next cycle, making the shock oscillation overstable. 

As the center of the primary vortex moves away from the symmetry axis during its last phase 
of the overstable oscillation {t/t 

cross ~ 465 — 492), it becomes less gravitationally bound and can 
thus be more easily influenced by the background flow. For all the supersonic models we run in this 
paper, we find that whenever the vortex rises above a half of the accretion radius, also known as the 
Hoyle-Lyttleton radius, defined by ta = 2GMp/Vp , it is swept away by the background flow in the 
negative s-direction (Fig. [6]). When the primary vortex is located outside 0.5rA, the counterstream 
associated with the vortex is weak and occurs in the far rear side of the perturber, unable to pass 
through the center of the perturber. The counterstream is instead resisted by the flow moving 
downstream right across the perturber and pushed up in the lateral direction (Fig. dig). This 
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develops another vortex with an opposite sense of rotation to the primary one, capable of pushing 
the latter away from the perturber (Fig. HFi). With the primary vortex carried away downstream, 
the flow near the perturber reaches a quasi-steady state in which shocked subsonic gas moves almost 
parallel to the symmetry axis. Since the associated kinetic energy is much smaller than the thermal 
and gravitational energies of the gas, the density distribution around the perturber becomes nearly 
hydrostatic, as will be shown below. 

4-1-2. Subsonic Cases 

Figure [7] displays density and velocity structures of a nonlinear subsonic model with ^ = 20 
and M = 0.5 in a conioving frame with, the perturber. Snapshots at ^/^cross — 

60, 150, and 600 

are shown. Sound waves launched from the perturber at t = propagate radially outward into the 
surrounding medium, forming a spherical causal region within which the medium is affected by the 
sonic perturbations. Since the perturber is spatially extended, however, the boundary of the casual 
region is not as sharp as in the case of a point mass, although the most dominant perturbations 
still come from the perturber center. Unlike supersonic cases, this model always involves subsonic 
flows and never produces a shock. Nevertheless, the overall flow pattern and late-time density 
structure near the perturber of this model is very similar to those in the postshock subsonic regions 
of supersonic models. First of all, the strong gravitational pull forms a counterstream and an 
associated vortex ring near the symmetry axis (Fig.[7]a). The counterstream moving in the upstream 
direction interacts with the incident flow (Fig. [7]6). The gas at the interface is pushed up toward 
the high-i? regions and then carried downstream, creating small vortices with low density (Fig.[7]c). 

The primary vortex slowly rises in the i2-direction due to buoyancy, and merge with the small 
vortices. Since this model does not contain a shock that would conflne the region of influence 
and since the causal region keeps expanding at a sonic speed, the flows in the high-i? regions are 
almost parallel to the symmetry axis. As a result, the primary vortex keeps rising as there is no 
momentum input in the background flow capable of pushing it back toward the symmetry axis. At 
the end of the run (t/tcross = 600), the primary vortex in this model arrives at (s, R) = (Sr^, 28rcj), 
corresponding to O.lSrA- It is uncertain whether the primary vortex will be carried downstream 
when it goes beyond r = 0.5rA in a manner similar to supersonic cases. At any event, the density 
distribution close to the perturber is well described by the condition of hydrostatic equilibrium at 
late time (see §4.21 below). 

4.2. Quasi-steady Density Wakes 

Figure [8] illustrates changes in the quasi-steady density wakes with varying A on the s-R plane 
for supersonic models with flxed M = 1.5 at t/tcmss = 600. The left panels show large-scale views 
of the wakes at —1400 < s/rg < 100 and < R/rg < 600, while the region near the perturber with 
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—60 < s/vg < 20 and < R/rg < 60 is enlarged in the right panels. In each panel, the perturber 
is located at s = i? = 0, and the black line connecting the points {s,R)/rs = (0,0) and (—60,54) 
marks the boundary of the Mach cone characteristics of the linear density wake for M = 1.5. When 
^ <C 1, the sonic perturbations are too weak to produce a shock, and the high-density ridge of 
the wake follows the Mach cone fairly well, although it is broadened due to the extended mass 
distribution of the perturber. As A increases to unity, sonic perturbations even outside the Mach 
cone attain substantial amplitudes enough to induce a bow shock that is attached to the center 
of the perturber (within the resolution limit). Since the density wake is effectively shifted toward 
the perturber compared with the linear case and still located preferentially at the rear side of the 
perturber, the resulting drag force will be larger than the linear counterpart (see ^4.4p . 

As A increases further, the shocked material gathered around the perturber begins to build up 
a strong pressure barrier which the incident flow cannot easily penetrate. This naturally makes the 
shock detached. Figure [8j*,c shows that the postshock density distribution around the perturber 
with A= 10 or 20 is almost spherically symmetric, indicating that the kinetic energy of the gas there 
is much smaller than the thermal and gravitational potential energies. To check if this is indeed 
the case, we plot in Figure [9] the density profiles along the symmetry axis and the -R-axis from the 
center of a massive perturber with A = 20; both supersonic and subsonic models with A4 = 1.5 
and 0.5 at t/tcross = 600 are presented. Also shown as dotted line is the density distribution under 
the assumption of hydrostatic equilibrium 

^ 1/(7-1) 

I , (7) 

where r = (s^ + i?^)^/^ denotes the distance from the perturber center and po and are the density 
and the adiabatic sound speed at r = 0, respectively. 

For both supersonic and subsonic models, the density distributions at r/r^ ^10 along the 
z^-, z~-, and i2-cuts are virtually identical to each other, and are in remarkable agreement with 
the predictions of equation ([TD with po/Poo = 40 and ao/aoo = 3.8 for M. = 1.5, and po/poo = 52 
and ao/aoo = 3.8 for M. = 0.5. The sharp drop offs of the density in the supersonic model near 
r/vs = 12 and 32 along the z~^- and -R-directions, respectively, are of course due to the bow shock. 
In the subsonic model, the presence of the primary vortex makes the local density decreased at 
R/vs = 28. Note that the quantity c? j p'^~^ measures the specific entropy and thus is conserved in 
an adiabatic flow without involving a shock, as is the case in the subsonic model. For the supersonic 
model, however, j p^~^ is increased from unity to 1.34 because of a shock jump. This corresponds 
to the shock Mach number of 2.3, which is larger than M. = 1.5 because the flow is accelerated 
by the perturber even before experiencing the shock. In the regions with r/r^ > 10, the density is 
overall larger along the z^- than z'^-directions, providing non-vanishing drag forces. Nonetheless, 
the presence of hydrostatic cores in the density wakes of massive supersonic perturbers makes the 
drag forces smaller than the linear results. 

While the density wake near a massive perturber with ^ > 1 is quite different from the linear 



(7 - ^)Aa 



2 

oo 



(r2 + r2)i/2 
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counterpart, the distant part of the wake is more or less the same. Figure [10] plots the exemplary 
profiles of the normalized perturbed density, ajA, along the cuts with i? = and Rjrg = 200 for 
models shown in Figure [8j Note that a/ A is nearly independent of A in regions far away from 
the perturber (e.g., s/rg < —900 region in Fig. [TOb and s/rg < —500 region in Fig. [TOb): in these 
regions, the gravitational potential perturbations are weaker by more than two orders of magnitudes 
than at the perturber location and thus locally in the linear regime. Even with low amplitudes of 
local perturbations, however, the region with —900 < s/rg < —100 along the symmetry axis behind 
the perturber has the nonlinear density wake that deviates considerably from the linear case. This 
is because the gas flowing in this region was already affected by strong gravitational potential in 
the upstream region, and has a diverging velocity field that reduces the perturbed density. Small 
fluctuations of nonlinear density wakes apparent in Figures [8] and [10] near the R = axis is thought 
of as arising from sonic perturbations induced by the primary vortex and its oscillations discussed 
in SU 



^=T7^' (8) 



4.3. Detached Shock Distance 

We have shown in the previous subsection that a sufficiently-massive supersonic perturber 
generates a density wake that is characterized by a bow shock standing ahead of the perturber and 
a surrounding hydrostatic envelope. Figure [5] shows that the quasi-steady value of the detached 
shock distance 6 is larger for models with larger A or smaller Ai{> 1). To quantify the dependences 
of 6 upon A and M, we introduce the nonlinearity parameter 

A 

and plot in Figure [11] the normalized shock distance against rj. The various symbols give the mean 
values of 5 temporally averaged over i/tcross > 50 that ignores the initial wake-development phase. 
The standard deviations of 6 are also indicated by errorbars. The numerical results are remarkably 
well described by the two simple power laws: S/vg = 2{ri/2)'^'^ for 0.7 ^ rj < 2 and 5/rg = rj for 
?7 > 2. 

The behavior of 5 with rj can be qualitatively understood as follows. For a very massive per- 
turber, the kinetic energy of the incident flow along the symmetry axis is almost entirely converted 
to the thermal energy of the postshock flow that supports a hydrostatic envelope against the grav- 
itational potential of the perturber. When the shock is strong {M ^ 1), the postshock thermal 
energy proportional to a^A4^ balances the gravitational potential energy —GMp/5 at the shock 
location, resulting in 5/rg oc A/ M."^. In the limit of — > 1, the flow does not in principle produce 
a shock, corresponding to (5 ^ oo. In practice, the gravitational acceleration is able to turn an 
incident nearly-transonic gas into a weakly supersonic flow, but the shock that may form is located 
very far away from the perturber anyway. In view of the detached shock distance, rj may be a bet- 
ter indicator of nonlinearity than A] for fixed A, the shape of a density wake due to a supersonic 
perturber becomes similar to the linear counterpart as M increases. 
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It is well known from hydrodynamic experiments and corresponding theories that a supersonic 
flow over a solid object with a blunt nose develops a detached bow shock when the nose angle 
is larger than the maximum angle allowed by the postshock flow (e.g., iLiepmann &: Roshkd 119571 . 
and references therein). The shocked gas becomes subsonic and slowly adjusts its velocity as it 
flows downstream to meet the boundary conditions at the surface of the object. In our simulations, 
the incoming flow toward a massive perturber recognizes the hydrostatic envelope as a spherical 
obstacle. Since the nose angle of a spherical body with respect to the incident flow is 90 degrees, 
the shock must be detached. 

Even though near-hydrostatic envelopes that form in our simulations are not entirely impen- 
etrable, we want to measure their effective sizes as perturbing obstacles. This can be achieved by 
comparing our numerical results for the detached shock distances with those from non-gravitating 
hydrodynamic theories (or lab experiments). Assuming that a bow shock ahead of a spherical body 
with radius Rs has a spherical shape near the symmetry axis. iGuyl (jl974l ) showed that the standoff 
distance of the shock is approximately given by 



where 



with fi'^ = ((7 



_5_ _ 

Rs 

1 _ 1 

K{M) ~ 2 

1)M^ + 2)/(2-yM '^ - 



4:{M^ - 1) 



+ 1 



l/i2K) 



1, 



1 + 



2 

7 + 1 fi 



2;u + l + 



(9) 



(10) 



7 + 1). Equation ([9l) has proven to explain the experimental 



data qu ite well ("e.g.. iHeberle et al.l Il950l : ISchwartz &: EckermanI Il956l : Ivan Dyke &: MiltonI Il958l : 
see also ISchreierlll982l ). Figure [T2] plots as a solid line 5/Rs from equation ([9]) with 7 = 5/3 as a 



function of M. Also plotted as various symbols are our numerical results for the ratio of 6 to the 
BHL radius, r^uL = GMp/{Vp + a^) for models with r] > 2, averaged over t/tcross > 50. Again, 
errorbars indicate the standard deviations. A rough agreement between the two results suggests 
that the BHL radius can be a useful measure of the effective size of the hydrostatic sphere. 



4.4. Gravitational Drag Force 

For a given wake distribution p{x, t) at time t of a perturber, it is straightforward to calculate 
the DF force on it by performing integration in equation ([5]). Figure [13] plots temporal changes of 
the DF forces normalized by 47rpoo(GMp/aoo)^ for models with differing A, but with fixed M. = 1.5, 
over the course of the wake evolution. The dotted line corresponding to the linear DF force (eq. [T]), 
with rjnin given in equation ([6]), closely follows the numerical results for A = 0.01, showing that the 
linear drag force increases logarithmically with time. The drag forces for nonlinear cases with high 
A also have a similar logarithmic time dependence, although they fluctuate for a while in response 
to the oscillations of primary vortices as well as detached bow shocks before a quasi-steady state is 
attained; the fluctuation amplitudes are typically ~ 4 — 16%, with a smaller value corresponding 
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to larger A and A4. Note that for M = 1.5, the normahzed drag force decreases with increasing 
A, indicating that the nonhnear effect makes the DF force smaller than the linear estimate. This 
is because a higher value of A implies a correspondingly larger detached shock distance, and a 
hydrostatic envelope with front-back symmetry near the perturber contributes a negligible amount 
to the net DF force. 

The dimensionless drag forces at t/tcmss = 600 when the wakes are in a quasi-steady state 
are given in Figure [T^ for various models with different A and Ai. For all the supersonic models, 
the DF force is a decreasing function of A. The reduction of the normalized DF force is larger 
for models with 7W ~ 1 than highly supersonic models. The Mach number corresponding to the 
maximum drag force shifts from unity to ~ 1.5 as the wake becomes highly nonlinear. For subsonic 
models, on the other hand, the nonlinear drag forces for ^ = 20 and 50 show some fluctuations 
(represented by errorbars) associated with slowly-evolving vortices present in the wakes, but their 
respective time-averaged values are very close to the drag forces in the linear regime. In fact, the 
similarity between the linear and nonlinear drag forces on subsonic perturbers is expec ted since 



the li near density wakes also possess front-back symmetry in the vicinity of the perturber (jOstrikei 



1993 ). Since the subsonic DF forces are dominated by the far field where perturbations are weak 



regardless of the perturber mass, the linear results should be valid even for very massive perturbers. 

The gravitational drag force certainly depends on both A and A4, but the discussions given 
above suggest that it may be through the nonlinearity parameter rj. To check this, we plot in 
Figure [l5] the ratio of the nonlinear DF force F to the linear prediction Fn^ as a function of rj. 
Again, various symbols give temporal averages of F/Fw^^ over t/tcross > 50, and their standard 
deviations are indicated by errorbars. For r] < 0.7 with which a bow shock that barely forms 
is attached to a perturber, F/F^^ ~ 1. When r] is increased to ~ 0.7-2, the shock becomes 
detached, but its standoff distance is not so large. In this case, most of the material in the wake 
is still located behind of, but closer to the perturber in comparison with the linear wake (see, e.g., 
Fig. [8]6), resulting in the DF force slightly larger (by less than 20%) than the linear counterpart. 
In highly nonlinear cases with rj > 2, however, the presence of a large hydrostatic envelope makes 
the drag force reduced considerably. For r] > 2, the numerical results are well fitted by 

-0.45 

,2. 



F = i^lin (?) . (11) 



Figure [T5] also plots the gravitational drag force from the hydrodynamic simulations of IShima et al 



(jl985l ) for BHL accretion flows, which are consistent with our numerical results. We defer to 



more detailed discussion of our results in connection with the BHL flows. 

To ascertain that the reduction of the DF force in highly nonlinear supersonic cases is really 
caused by the presence of spherically-symmetric hydrostatic envelopes near the perturbers, we 
calculate the drag force by imposing a cut-off radius r^[^ such that only the regions in the wake 
with r > Tmin participate in the force evaluation (eq. [5])- Figure [16] plots the resulting dependence 
of F upon rinin for the model with ^ = 20 and M = 1.5 at t/tcmss = 600. The vertical dashed 
line marks the location of the bow shock along the symmetry axis, while the dotted line indicates 
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a slope of — 1. Note that the drag force is independent of rmin for Tmin < S, clearly demonstrating 
that the hydrostatic sphere surrounding the perturber has a negligible contribution to the net drag 
force. Figure [TBI shows F oc — ln(rmm/''"s) for large r^in, analogous to the linear cases (se e eq. [1]). 
From the study of BHL accretion flows to large gravitating bodies, IShankar et al.l (|l993l ) similarly 
found that the drag force declines logarithmically as the size of the accretor increases. 



4.5. Resolution Dependency 

Finally, we remark on the effect of numerical resolution on the DF force. Figure [T71 plots the 
time evolution of the detached shock distance as well as the drag force on a perturber with ^ = 10 
and A4 = 1.5 from models with different resolution: Az/rg = 0.1, 0.2, 0.4, and 0.8, where Az is the 
grid spacing. In models with Az/vg = 0.4 and 0.8, the shock has a larger standoff distance without 
noticeable oscillations and the DF force is correspondingly smaller than in models with higher 
resolution. Compared with the Az/r^ = 0.2 model, the model with Az/r<j = 0.1 presents the shock 
oscillations with higher amplitudes and arrives at quasi-s teady equilibrium ear lier. This resolution 



dependence of the shock oscillations was also noticed by iMatsuda et al.l (|l989l ) for adiabatic BHL 
flows onto an absorbing perturber. The resolution study shown in Figure [17] suggests that our 
numerical results for the DF forces are reliable as long as 5 or more grids per are taken. 



5. SUMMARY AND DISCUSSION 

DF of bodies orbiting in a gaseous medium is of great importance in various astronomical 
systems ranging from protoplanetary disks to galaxy clusters. In previous analytic studies of DF, 
it has been assumed that the mass of a moving object is small enough for the induced density wake 
to have low amplitudes and be thus in the linear regime. However, there are many astronomical 
situations such as in a merger of black holes near a galaxy center and migration of protoplanets, 
where a perturber is so massive that the induced wakes are well in the nonlinear regime. In this 
paper, we use numerical hydrodynamic simulations to explore nonlinear gravitational responses 
of the gas to, and the resulting drag force on, a perturber with mass Mp moving straight at 
velocity Vp in an initially static, uniform background with density poo and adiabatic sound speed 
Ooo. The perturber is represented by a Plummer sphere with softening radius r^. Unlike in the 
BHL problems, the perturber in our models does not possess a defined surface through which gas 
is accreted or reflected. Assuming an axial symmetry, we solve the basic equations on the (i?, z) 
plane, where R and z denote the direction perpendicular and parallel to the motion of the perturber, 
respectively. Our numerical models are completely characterized by two dimensionless parameters: 
A = GMp/ {a^Ts) and M. = Vp/aoo- To study DF in various situations, we run as many as 58 
models that differ in A and A4. Our standard models have 5 zones per r^, but we also run models 
with different resolutions to ensure that the density and velocity wakes are fully resolved. 
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For supersonic models {M > 1), we find that a massive perturber with ^ ^ 1 produces 
a bow shock ahead of it through which the incident supersonic flow becomes subsonic. In the 
beginning, the postshock flow in the nonlinear cases develops transient features such as a primary 
vortex ring and an associated counterstream near the symmetry axis, which causes the shock to 
oscillate around its equilibrium position. The shock oscillation added by the counterstream is 
overstable, amplifying the amplitudes of its oscillation and the vortex movement in the {R, z)- 
plane (see §4.ip . The vortex ring is eventually shed downstream from the perturber, leaving the 
nearly-stationary bow shock and a quasi- hydrostatic envelope that surrounds the perturber. On 
the other hand, subsonic models with 1 without involving a shock retain a primary vortex 

and many other small-scale structures until the end of runs. Nevertheless, strong gravity makes 
the density distribution near the perturber still very close to that under hydrostatic equilibrium. 

By comparing the numerical results from various supersonic models with differing A and A^, 
we find that the simulation outcomes such as the detached shock distance and the gravitational DF 
force are very well quantified by a single parameter defined in equation ([8]). When rj < 0.7, the 
system is in the linear regime where a Mach cone is bounded by low-ainplitud e Mach waves, and the 
drag force is just the same as the analytic linear value F\i^ of lOstrikeii (|l999l ). When rj is moderate 
at 0.7 ^ rj < 2, the Mach waves turn into a bow shock that is weakly detached, with the standoff 
distance 6 ~ 2{rj/2)'^'^rs from the perturber. In this case, the density wake is slightly shifted toward 
the perturber compared with the linear counterpart, causing the drag force to be larger than F^^ 
by a factor of less than 1.2. In the highly nonlinear regime with r] > 2, however, the detached shock 
distance behaves as 6 = r]rs, and the nonlinear drag force F is given by F/Fhyi = {r]/2)^^''^^ . The 
reduction of the drag force compared with the linear value is because the density wake close to the 
perturber is in near-hydrostatic equilibrium and thus contributes a negligible amount to the net DF 
force. Since front-back symmetry notable for nonlinear wakes exists also in the linear wakes with 
M < 1, the nonlinear drag force on a massive subsonic perturber is similar to the linear prediction. 

As mentioned in 311 our model simulations differ from those of the BHL accretion flows in terms 
of the boundary conditions. In our models, a perturber simply provides a gravitational potential for 
the background gas and does not hold any surface, while models for the BHL accretion considered 
a defined surface through which the gas is accreted or refiected. The different boundary conditions 
might lead to different evolution and structure of wakes. For instance, the BHL accretion fiows 
around a perturber present nonlinear features including a detached bow shock, vortex oscillations, 
and counterstreams in t he forward /backwa, r d directions, just a s in our models, if the perturber is 
not perfectly absorbing (jShima et al.l Il985l : iFrvxell et al.l 119871 : iMatsuda et al.l Il989l ) . When the 
perturber has a totally absorbing surface, on the other hand, the wakes are relatively quiescent, 
fiowing nearl y spherically into the perturber that a,bsorbs the angular momentum , carried by the 
accreting gas (jShima et al.lll985l : lFrvxell et al.lll987l : lKoide et al.lll99ll : lRuffertlll994l ). In the case of 
an extremely small perturber compared with the accretion radius, even an absorbirig bod y produces 



violent features since it is not able to accept the whole accreting gas (IKoide et al 



1991 



Despite these differences, the detached distances of bow shocks and the resulting DF forces 
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appear to be not so sensitive to the boundary conditions adopted. Figure [18] plots the detached 



shock distances inferred f r om the results of two-dimensional axisym metric simulations (jShima et al 



19851 : iFrvxell et al.1 Il987l : iMatsuda et alJ Il989l: iKoide et alJ Il99lh and the tabulated values from 
three-dimensional nonaxisymmetric runs (jRuffertl 11994 ) for the BHL accretion flows, by taking 
equal to the radius of the perturber. The BHL results under the reflecting boundary conditions 
are larger than our results by a factor of only ~ 1.2, while those under the absorbing boundary 
conditions are smaller by a factor of ~ 0.4. Since the work on the BHL accretion mostly focused 
on mass accretion rates, only a few of them evaluated the gravitational drag forces. In Figure 1151 
we plo t using star symbols the drag forces from the adiabatic runs with 7 = 5/3 in IShima et al 



(jl985l ). which are the only data we acquired from the literature that allow a reliable comparison 
with our result^. Note that the drag force is larger on a purely-abso r bing perturber for which the 
detached shock distance is smaller. The DF forces from IShima et al.l (|l985l ) are roughly consistent 
with, and follow a similar trend to, our results. This is probably because even though density wakes 
are nowhere close to being hydrostatic in BHL accretion flows, they somehow maintain front-back 
symmetry that makes the net DF force smaller. 

While we have considered in this paper only an adiabatic gas with index 7 = 5/3, we note 
that the DF force may depends sensitively on 7. Simulations of the BHL accretion flows reported 
that as 7 decreases, the detached shock distance decreases and th e dens i ty wake becomes shaped 



increasingly i nto a cone-like structure similar to the linear wake (jHunt 



Ruffert 



1995 



197S 



Shima et al 



1985 



19961 ). This is presumably because bow shocks in lower-7 models should be stronger 



in order to compensate for the diminished pressure in supporting the postshock ga s against the 



pertu rber's gravity. Consequently, drag forces are larger for the gas with smaller 7 (IShima et al, 
1983). 



Escala et al.l (j2005l ) carried out numerical simulations of a DF-induced merger of supermassive 
binary black holes due to an adiabatic gas with 7 = 5/3. Their Figure 9 presents the effect of 
the black hole mass on the evolution of the binary separation. From this figure, we infer that the 
orbital decay time scales with the black hole mass as Tdecay oc M~^'^ in their simulations, which is 



somewhat shallower than the linear expectation of T^ecay oc Mp^ (e.g.,l0strikeiJll999|). On the other 
hand, our results predict Tdecay oc M~^-^^ if the perturb er is sufficient l v ma ssive, suggesting that 
the delayed orbital decay of supermassive black holes in lEscala et al.l (120051 ) is partly due to the 
nonlinear effect. Of course, there are many other factors that may change Tdecay While we consider 
quite an ideal situa t ion in which a perturber is moving straight in a static, uniform medium, the 
gas in lEscala et al.l (j2005l ) is distributed in a rotating, radially-stratified, self-gravitating disk and 



the black holes follow curvilinear, possibly eccentric, orbits. In addition, the Mach number of the 
black holes in their models is likely to vary during the decay toward the orbital center, which may 
also modify the decay time. Since these compact objects have a very small size, it is also an issue 



Although IShankar et al. I l|l993t ) and lRuffertI ([199J) also gave the drag coefficients the perturbers are too large in 
the former and the units are uncertain in the latter to be compared with our results. 
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whether numerical models resolve detached bow shocks, which is crucial in evaluating the drag 
force accurately. We will discuss these in subsequent work. 
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A. Linear Solution for the Velocity Wake 

Usi ng the Fourier transf orm for the space variables and the Laplace transform for the time 



variable, lJust Sz Kegell (|l99d ) derived an analytic solution for the velocity field in the linear wake, 
but unfortunately the resulting expression (their equation [47]) contains an typographical mistake. 
Here we rederive the expression using Fourier transforms for both space and time variables. We 
consider a point-mass perturber with mass Mp on a straight-line trajectory through a uniform 
medium with density poo and adiabatic sound speed aoo- The basic equations of hydrodynamics 
for an inviscid, nonmagnetic, non-self-gravitating gas are 

1^ + V . {pv) = 0, (Al) 

and 

dv 1 

+ Vt; = — VP- V^-ext, (A2) 
ot p 

where $ext is the gravitational potential of the perturber with density pext, satisfying the Poisson 
equation V^^>ext = 47rG/9ext- Assuming that the density and velocity fields induced by the perturber 
have small amplitudes, we linearize equations (|Aip and (|A2p to obtain 



dt 

and 



^ + V • ^ = 0, (A3) 



— + aLVa = -V^ext, (A4) 



dt 

where a = p/poo — 1 is the dimensionless perturbed density. 

Taking Fourier transformation of equations (|A3p and (IA4p , we have 



,2 ^'^\^ 47rG_ 

k - — " = Pext (A5) 



at 



and 

V = -p-") (A6) 
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where symbols with and without hat denote a Fourier pair defined by 



(27r)^ 




and 



^{k,uj) = JJJJ ^{x, f)e-^(^-^-'^*) d^x dt, 
with ^ representing any physical quantities. 



(A7) 
(A8) 



For a perturber introduced at t = and moving at a constant speed Vp along the positive 
z-direction, pext = Mp 6{x — Vpte^) 'H{t), where 6{x) is the Kronecker delta and TC{t) is a Heaviside 
step function. This results in 



Pcxt 



2Tr5{uj - k,Vp) + — — 



(A9) 



The second term in e quation (IA9D is originated from the step function, which vanishes if a steady 



state is assumed (e.g.. iFurlanetto &: Loebll2002l ). 



Substituting equation ()A9P into equation (jASp . and taking an inverse Fourier transform of 
the resulting expression for S?, one can obtain the solution of the perturbed density for finite-time 
pe rturbations. The de tailed procedure is similar to that of the Fourier-Laplace transform method 
of I Just &: Kegell (|l99Cll ). so that we simply write the result 



GMr, 



a 



Ti 



(s2 + i?2(i__yv^2))i/2' 



where s 



-Vpt and Ti 



2, 1, and in region I {R^+z^ > a^t^, W < 



-sz. 



(AlO) 



s^ + R\l-M^) > 0, 



Ai > 1), region II (i?^ + z'^ < a^t'^), and any other region, respectively. Note that equation (jAlOp 
is identical to equation (10) of [OstrikeiJ (|l999l ) based on the retarded Green's function technique. 



Combining equations (1A5 
transform, we have 



A6h , and (lA9h together in favor of v and taking its inverse Fourier 



GMp 




2ir5{uj — kzVp 



+ 



i(k-x- 

k^ 



^^U^kduj. (All) 



The second term in the integrand of equation (lAlip has three simple poles on the real uj axis 
at uji = kzVp and a;2,3 = ±kaoo- To ensure the proper analyticity, we must displace them by 
an infinitesimal amount into the upper or lower half w-plane. Since our problem requires v = 
everywhere for t < 0, the proper choice of the pole displacements should be oji into the upper half 
w-plane and 102,3 to the lower half cj-plane. For t < 0, then, a contour of integration consisting of 
the real to axis and an infinitely-large semicircle in the upper half w-plane encloses only wi, and 
the calculus of residues guarantees that after the integration over iv, the first and second terms in 
the integrand of equation (|Alip become identical with each other with opposite sign, resulting in 
no velocity perturbations before the introduction of the perturber. 
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For t > 0, the contour that goes along an infinitely-large semicircle in the lower half cj-plane 
instead contains the poles at 012,3 • The remaining steps in the integration with respect to uj is quite 
straightforward to yield 



GMp 



k + k.M 



-ik^Vvt 



k — kzM 



k^ 



X j3 



(A12) 



To carry out the integration in equation ()A12p . it is convenient to use spherical coordinates in the 
fc-space as 

k = k {s'mO cos{{p + Q , sin 9 sm{ip + (),cos9), 
which is oriented relative to the space coordinates 



X = {R cos C, R sin C, z) 
satisfying k ■ x = k (z cos 9 + R sin 9 cos 99). We then obtain 



(A14) 



V 



GMp 



^ik(z cos 8+aaat) ^iks COS 8 

l + Mcos9 



^ik{z cos d—aoot) ^ikscost 

I- Mcos9 



^ikR sin 6 cos ip 



X (sin^ 9 cos{(p + C)ex + sin^ 9 sin{(p + C)ey + sin 9 cos ^e^) dk d9 dip, (A15) 



where ex, Cy, and are the unit vectors in the x-, y-, and z-directions, respectively. Without 
loss of generality, we take C = 0, in which case ex and ey correspond to the unit vectors in the 
radial (eR,) and azimuthal (e^) directions, respectively. It then follows that the azimuthal velocity 
vanishes (i.e., Vip = 0) since this component in the integrand is an odd function of (p. 

For the velocity component in the direction of motion, we substitute 9 = tt — 9' to the second 
term inside the square bracket of equation ()A15j) to obtain 



GM„ 



vra 



oo JO 



d9 / dk 



sin 9 cos 9jQ{kR sin 9) 







1 + Mcos9 

Using [G6.671.8jl, inserting the substitutions 

{R^ + z^) cos 9 + zaoot 
R{R'^ + z2 - a^t2)i/2 



{cos[k{z cos 9 + aooi)] — cos(A;s cos 9)} 



smry 



and 



(i?2 + s2)l/2 

sm 7] = cos ( 

R 



(A16) 



(A17) 



(A18) 



for the first and second terms of equation ()A16p . respectively, and applying [G2.551.3] for the 
r/-integrations, we obtain 



GMp 



(ii2 + s2)l/2 (s2+^2(l__/^^2))l/2 (^2 + ^2)1/2 



(A19) 



■^[G- • •] refers to the number of formula in the integral table of 



Gradshtevn &: Rvzhik 



mm . 
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where Ti =2, 1, and T2 = 1, 0, 1 for region I {R^ + z'^ > a^t^ < -sz, + R'^{1 - M'^) > 0, 
A4 > 1), region II (i?^ + < a^t^), and any other region, respectively. 



For the velocity component in the radial direction, the use of [G6. 671.1] similarly leads to 



VR 



GMp 



zT, 



(i?2 + s2)l/2 (52 + ^2(l__yV^2))l/2 (i?2 + ^2)1/2 



(A20) 



The results of our time-dependent linear analyses for density and velocity fields are confirmed by 
direct numerical simulations for a low-mass perturber with A = 0.01 presented in ^ 

In terms of our notatioE0, equation (47) of Just Sz Kegel ( 1990 ) can be rewritten as 



GMp 



is,-R) 



{s,-R)Ti 



{z,-R) 

(i?2 + s2)l/2 (s2 + i?2(l__y^2))l/2 (^2 + ^2)1/2 



(A21) 



Comp arison of equation (|A2ip with equations (|A19p and (jA20p shows that equation (47) of lJust &: Kegel 
(jl990l ) takes T2 = 1 everywhere, which is clearly incorrect. 
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Fig. 1. — Color-scale distributions of the perturbed density a (top), parallel velocity Vz (middle), 
and perpendicular velocity vr (bottom) to the line of motion of an extended perturber with A = 0.01 
and M = 1.5. Colorbars label a, Vz/aoo, and VR/aoo from top to bottom. The black contours 
represent the results of the time-dependent linear perturbation theory for the corresponding point- 
mass perturber. The perturber initially located at (z, R) = (0, 0) has moved to (Maoot, 0) at time 
t along the z-axis. 
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Fig. 2. — Distributions of the perturbed density along the cuts at R/aoot = 0.20 {top) and z/aoot = 
0.92 (bottom) indicated as dotted Hues in Fig. [1] of a model with A = 0.01 and A4 = 1.5. Open 
circles representing the simulation results deviate considerably from the results a (dashed line) of 
the linear perturbation theory for a point mass, but are in excellent agreement with Oext (solid 
line) obtained from the convolution of a with the extended mass distribution. 
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Fig. 3. — Dimensionless DF force on a low-mass perturber with A = 0.01 against the Mach 
number A4 at t/tcross = 300. The sohd hne is the best fit of the hnear force formula (eq. [T]) with 

'"mill — 

our simulation results {open squares). 
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Fig. 4. — Temporal evolution of the density structure (logarithmic color-scale) and velocity field 
(arrows) generated by an extended perturber with A = 20 and A4 = 1.5 in the frame where the 
perturber located at (s = z — Vpt, R) = (0, 0) is stationary. A black semicircle marks the softening 
radius of the perturber, while red curves draw the sonic lines where the gas speed is equal to the 
local sound speed. Colorbar labels \og{p/ poo — !)• 
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Fig. 5. — Time evolution of the detached shock distance 5 measured along the symmetry axis from 
the perturber center for various models. For each model with A < 50, 5 initially increases rapidly 
with time, experiences quasi-periodic oscillations, and then saturates to a constant value. The time 
to reach a quasi-steady state and the value of 5 at saturation become smaller as A decreases or M. 
increases; the model with A = 100 and M. = 1.5 does not attain a quasi-steady state until the end 
of run. 
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Fig. 6. — Trajectory of the center of the primary vortex on the s-R plane for a model with A = 20 
and A4 = 1.5. Filled circles marked by alphabets a-h indicate the vortex locations at time epochs 
corresponding to the snapshots shown in Fig. HI The dotted curve draws the region where the 
distance from the perturber equals a half of accretion radius ta = 2GMp/Vp . 
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Fig. 7. — Snapshots of the perturbed gas density (logarithmic color-scale) overlaid with the velocity 
structure (arrows) of a nonlinear subsonic model with A = 20 and A4 = 0.5 in a comoving frame 
with the perturber. Colorbar labels log{p/poo — !)• 
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Fig. 8. — Distributions of quasi-steady density wakes with varying A on the s-R plane at t/tcross = 
600, with the perturber placed at (s, R) = (0, 0). Left panels show large-scale views of a, while the 
right panels focus on a small section with —60 < s/rg < 20 and < R/rg < 60 near the perturber. 
All the models have the same M = 1.5. The black line in each panel traces the Mach cone formed 
by a low-mass point-mass perturber with the same Mach number. The perturber gathers more gas 
toward it as A increases. Colorbars label log{p/poo — !)• 
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Fig. 9. — Density profiles along the positive- 2; {solid), negative-z (dashed), and R (dot-dashed) 
axes from the perturber center for (a) a supersonic model with A = 20 and A4 = 1.5 and (6) a 
subsonic model with A = 20 and A4 = 0.5 at t/tcross = 600. In each panel, the dotted line gives the 
respective density distribution under the assumption of hydrostatic equilibrium, which is in good 
agreement with the numerical results for r/r^j < 10. 
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Fig. 10. — Distributions of normalized perturbed density aj A as functions of s along the (a) R = 
and (6) R/rg = 200 cuts for models shown in Figure El While a/ A near the perturber with A>1 
deviates significantly from that of the linear case, it becomes nearly independent of A in the regions 
far away from the perturber. For A> 1, a shows some fluctuations and becomes negative in some 
regions near the symmetry axis {R = 0) behind the perturber, since the flow in these regions already 
experienced strong perturbations in the upstream region. 
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Fig. 11. — Detached shock distance 6 as a function of the nonlincarity parameter rj = A/{M'^ — 1). 
Various symbols and errorbars give the means and standard deviations of 6 over time for t /tcross > 
50. Dashed Unes correspond to the broken power laws S/rg = 2(r//2)^-® for 0.7 < ry < 2 and S/vg = rj 
for 77 > 2 that describe the numerical data quite well. 
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Fig. 12. — Dimensionless detached shock distance as a function of the Mach number. Various 
symbols and errorbars give the means and standard deviations of 6 /tbrl 

over tiniG for ^/^cross 

> 50, 

where tehl is the BHL radius. The sohd curve plots equation ([9]) for (5/i?s , the ratio of the standoff 
distance of the shock to the radius of a spherical body in the non-gravitating hydrodynamic theory. 
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Fig. 13. — Temporal changes of the dimensionless DF forces for models with = 1.5. The dotted 
line is the prediction of equation ([T|) with rmin = O.SBJ^^'^rg, in good agreement with the numerical 
result for a low-mass perturber with A = 0.01. The drag forces increase logarithmically with time, 
although models with large A exhibit early fluctuations in accordance with shock oscillations. 
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Fig. 14. — Dimensionless DF forces for various models with different A. and A4 at t/tcmss — 

600. 

The sohd line draws the linear drag force (eq. [1]) with rmin = O-SbAd^'^rg. For supersonic models, 
the drag force decreases with increasing A. Subsonic models do not reach a quasi-steady state and 
show some fluctuations, but their mean drag forces are similar to the linear estimates with the 
same A4. 
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Fig. 15. — Ratio of the nonlinear to linear DF forces for supersonic models as a function of the 
nonlinearity parameter r] = A/{Ai'^ ~ !)• Various symbols and their errorbars indicate the means 
and standard deviations of F/F\\^ over time for t/t^j-oss 

> 50. When r/ < 0. 7 F/Fn^ ^ 1 while 
P/Fiin ~ {r]/2)~^-^^ for rj > 2. Also plotted as star symbols are the results of IShima et al.l (jl985l ) 
for BHL flows onto accreting or reflecting bodies. 
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Fig. 16. — Dimensionless DF force as a function of the cut-off radius rmin, interior of which is 
excluded in the force evaluation, for a model with A = 20 and M. = 1.5 at t/tcross = 600. The 
vertical dashed line marks the detached shock distance oi 5 = 12rs, while the dotted line indicates 
a slope of —1. The drag force is nearly constant for rmin ^ <^ and decreases almost logarithmically 
at large rmin- 
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Fig. 17. — Time evolution of (a) the detached shock distance S and (b) the dimensionless DF force 
F for ^ = 10 and A4 = 1.5. The simulation results at four different resolutions are compared. 
Note that 6 and F are fully resolved if Az/vg ^ 0.2, where Az is the grid spacing. 
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Fig. 18. — Detached shock distances against rj from the published work on the BHL accretion flows. 
Open symbols denote the results when the matter is allowed to be absorbed to the perturber, while 
those under the reflection boundary conditions are given by filled symbols. For comparison, the 
broken power laws that fit our numerical results well are also shown as dashed lines. 



